Method and system for characterizing tumors

ABSTRACT

Embodiments of the present invention are directed to determining the spatial extent, aggressiveness, and other characteristics of various types of tumors, including glioma tumors that occur in brain tissue. Various embodiments of the present invention use parameterized computational models to characterize tumor growth and employ medical imaging technologies to generate images and other types of data from which values of parameters of the computational models are derived. Having obtained the parameters for a particular tumor, the extent of the tumor is estimated, with high accuracy, and other characteristics of the tumor are derived from the parameterized computational models.

TECHNICAL FIELD

The present invention is related to characterizing tumors from dataprovided by various medical technologies, including imagingtechnologies.

BACKGROUND OF THE INVENTION

There are many different types of tumors and cancers. Gliomas are themost common type of primary brain tumors and are thought to arise fromglial cells, or precursors of glial cells, that surround, support, andcooperatively interact with neurons in the brain. Because gliomasgenerally have already extensively invaded brain tissue before a patientis aware of the disease, gliomas are generally not curable.Glioblastomas are the most malignant and most commonly occurring type ofgliomas in adults, representing about 50 percent of all gliomas.Glioblastomas are distinguished by necrosis, or tissue death, within thecentral portion of the tumor or irregularly spaced between vasculartissues, and by a surrounding shell of tissue diffusely invaded byperipheral tumor cells and characterized by edema. The aggressivebehavior of glioblastomas is reflected in the nearly 100 percentfatality rate of patients suffering from this type of tumor withinapproximately two years following diagnosis, even after extensivemedical intervention, including surgery, radiotherapy, and chemotherapy.

Despite continual advancements in imaging technologies, glioma cellsgenerally invade far beyond the regions recognized as being abnormalusing various types of clinical imaging technologies, includingcomputer-aided tomography (“CT”), magnetic-resonance imaging (“MRI”),and positron emission tomography (“PET”). The extent of glioma-cellinvasion is generally greater than that assumed for currentradiotherapy-treatment planning. Currently, an additional shell oftissue invaded by glioma cells, having a thickness of approximately twocentimeters, is assumed to surround the volume of the tumor seen in CT,MRI, or PET images. Because the extent of glioma-cell invasion iscurrently inaccurately modeled, therapies based on understanding theextent of glioma-cell invasion are often misapplied, leaving certain ofthe glioma cells that have diffused away from the tumor mass untreatedor inadequately treated. The medical community continues to seekimproved methodologies for characterizing both the extent ofglioblastomas and other tumors as well as determining the aggressivenessof tumors and additionally characterizing the tumors in order to guidemedical interventions as well as provide accurate prognoses forpatients.

SUMMARY OF THE INVENTION

Embodiments of the present invention are directed to determining thespatial extent, aggressiveness, and other characteristics of varioustypes of tumors, including glioma tumors that occur in brain tissue.Various embodiments of the present invention use parameterizedcomputational models to characterize tumor growth and employ medicalimaging technologies to generate images and other types of data fromwhich values of parameters of the computational models are derived.Having obtained the parameters for a particular tumor, the extent of thetumor is estimated, with high accuracy, and other characteristics of thetumor are derived from the parameterized computational models.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1A-D illustrate a simplistic computational model for tumor growth.

FIGS. 2A-C illustrate aspects of a more general computational model fortumor growth, employed by embodiments of the present invention.

FIGS. 3A-D illustrate growth of a glioblastoma, each figure representinga two-dimensional cross-sectional slice through the tumor mass.

FIGS. 4A-D illustrate a cell-proliferation component of tumor growth.

FIGS. 5A-D illustrate diffusive migration of tumor cells within atissue.

FIGS. 6A-F illustrate tumor growth, in cross-sections through a tumorvolume, as monitored by two different types of medical imaging.

FIGS. 7A-D illustrate preparation of a lookup table from whichcorrections to computed

$\frac{\rho}{D}$values are obtained according to one embodiment of the presentinvention.

FIG. 8 illustrates determination of a scalar or tensordiffusion-coefficient field according to one embodiment of the presentinvention.

FIGS. 9A-E illustrate one method by which the scalar diffusion fieldD(r) is determined.

FIG. 10 illustrates an alternative method for estimating the parametersρ and D according to one embodiment of the present invention.

FIG. 11 provides a control-flow diagram for a routine “tumor-growthanalysis” that represents one embodiment of the present invention.

FIG. 12 provides a control-flow diagram for a routine “use model” thatrepresents one embodiment of the present invention.

FIG. 13 provides a control-flow diagram for the alternative,single-image method for determining ρ and D parameters for a tumor.

DETAILED DESCRIPTION OF THE INVENTION

Various embodiments of the present invention are directed tocomputationally modeling tumor growth and then using the computationalmodel, along with a small number of medical images of the tumor, tocharacterize the spatial extent of the tumor as well as other tumorparameters, including the aggressiveness of the tumor. The methods ofthe present invention are necessarily carried out by electroniccomputing systems controlled by software programs, because handcalculation would be far too slow, inaccurate, and error prone to beuseful for medical diagnosis and medical therapy planning. In fact, handcalculation would be practically infeasible, due to the enormous numberof calculations involved in modeling growth of a particular tumor. Aswith any computational-modeling undertaking, there are many differentapproaches to modeling tumors, with the useful approaches characterizedby computational feasibility and simplicity while nonetheless providingdesired levels of accuracy and reproducibility.

FIGS. 1A-D illustrate a simplistic computational model for tumor growth.In FIG. 1A, a small, nascent tumor 102 is located within a larger volumeof tissue 104, with the location of the nascent tumor 102 located at theorigin of a three-dimensional Cartesian coordinate system. The tumor ismodeled as growing in a spherically symmetric manner, illustrated inFIGS. 1A-C, with the radius of the tumor increasing over time. A plot ofthe tumor radius versus time, as shown in FIG. 1D, would generallyprovide a curve 110 to which a computational model is fit to generatethe radius as a function of time: r=ƒ(t). Unfortunately, a simplisticcomputational model, such as that described above with reference toFIGS. 1A-D, would almost certainly not provide sufficient accuracy toinform medical decisions. Tumors generally do not grow in a sphericallysymmetric fashion. Moreover, tumors in different patients, or atdifferent locations within a patient, may grow at different rates.Tumors are also frequently not located at a central position within aspherical volume of tissue, and may be constrained by the geometry ofthe tissue within which they grow. For all of these reasons, and formany additional reasons, a more complex computational model for tumorgrowth is needed for medical applications.

FIGS. 2A-C illustrate aspects of a more general computational model fortumor growth, employed by various embodiments of the present invention.First, the tissue within which the tumor grows may have an arbitrarygeometry, such as the ellipsoid 202 volume of tissue shown in FIG. 2A.Although an ellipsoid has a high degree of symmetry, including an n-foldand two 2-fold symmetry axes coincident with the x, y, and z axes inFIG. 2A, arbitrary volumes, including low-symmetry volumes, can beaccommodated by the computational models used in various embodiments ofthe present invention. Second, the nascent tumor 204 may be locatedarbitrarily within the volume of tissue, according to computationalmodels employed in various embodiments of the present invention. Asshown in the series of FIGS. 2A-C, the tumor may grow in anon-spherically symmetric fashion.

FIGS. 3A-D illustrate growth of a glioblastoma, each figure representinga two-dimensional cross-sectional slice through the tumor mass.Initially, as shown in FIG. 3A, the nascent tumor is small and roughlyspherical 302. However, as shown in FIGS. 3B-D, as the tumor grows, bothan internal necrosis-containing region 306 and an outer,edema-associated region 308 have an increasingly complex andlow-symmetry shape. These complex, low-symmetry tumor volumes are aproduct of the inherent anisotropy biological tissues. Often, forexample, tumor cells may migrate preferentially along vascular tissue orother internal boundaries, which may have complex geometries. Many othertypes of anisotropy occur in biological tissues. In brain tissue, forexample, the distribution of white matter and gray matter isanisotropic, and the rate of glioma-cell migration may differ in the twodifferent types of tissues.

A number of different processes contribute to tumor growth and spatialexpansion. One process driving tumor growth is the proliferation oftumor cells, which generally divide, or proliferate, at rates that farexceed the normal proliferation rates of non-cancerous cells. FIGS. 4A-Dillustrate a cell-proliferation component of tumor growth. FIG. 4A showsa cross-section of a small volume of tissue containing mostly normalcells, such as normal cell 402, represented by unfilled disks, and asingle malignant cell 404, represented by a filled disk. Over time, asillustrated in the sequence of FIGS. 4A-C, the malignant cell divides tocreate progeny cells, and progeny cells, in turn, divide to produce afinal, relatively high concentration of malignant cells within thevolume of tissue, as shown in FIG. 4C. Often, there is a maximumconcentration of malignant cells, k, at which point furthermalignant-cell proliferation ceases and necrosis ensues. A simpledifferential-equation model for the cell proliferation illustrated FIGS.4A-C is:

$\frac{\partial c}{\partial t} = {\rho\;{c\left( {1 - \frac{c}{k}} \right)}}$where

-   -   c=the concentration of malignant cells at points in a tissue        volume, and is therefore a scalar function over a bounded        Cartesian space (x, y, z) and time;    -   t=time;    -   k=maximum malignant-cell concentration; and    -   ρ=a net rate of malignant-cell proliferation.        Solution of this differential equation produces a function c(x,        y, z, t) exponential in t. FIG. 4D plots a solution for this        differential equation as a concentration of malignant cells, c,        with respect to time, t. The concentration initially increases        exponentially and then asymptotically approaches the maximum        concentration in the malignant cells k.

Another process that contributes to tumor growth is tumor-cellmigration, which can be represented as a diffusion process. FIGS. 5A-Dillustrate diffusive migration of tumor cells within a tissue. FIGS.5A-C illustrate a small cross-section of a small volume of tissuecontaining tumor cells. In each of these figures, two positions withinthe volume are represented by position vectors r₁ 502 and r₂ 504. InFIG. 5A, the tumor cells, represented as unfilled disks, such as tumorcell 506, are clumped together in a small, central portion of thecross-section of the small volume of tissue. Over time, as representedby the series of FIGS. 5A-C, the cells migrate randomly, ending up, inFIG. 5C, relatively widely dispersed throughout the small volume oftissue. Note that, in the diffusive migration process, the concentrationof malignant cells, c, decreases over time at position r₁, but increasesover time at position r₂. In the computational modeling of diffusionprocesses, the direction of a gradient vector computed at a point (x, y,z), or {right arrow over (r)}, determines, for a next short interval oftime, whether the concentration increases or decreases, due todiffusion, and the magnitude of the gradient vector is related to fluxof diffusing entities, for a next short interval of time. For example,FIG. 5D shows a topographic-map-like plot of concentration within atwo-dimensional cross-section of a tissue, where the concentrations areplotted using contour lines of constant concentration. Thus, FIG. 5Dshows a three-dimensional concentration surface above the x,y plane,where the concentration of malignant cells c=ƒ(x,y). The differentialoperator ∇ can be applied to a scalar function to produce a vectorfield, generating a gradient vector emanating from each point upon whichthe scalar function is defined. In the plot shown in FIG. 5D, forexample, the gradient vector (g=∇c) 510 at position {right arrow over(r)}, points in the direction of the steepest concentration decrease.

To model the diffusive migration process, the concentration of malignantcells is assumed to be a scalar function of location and time:

$\begin{matrix}{c = {f\left( {x,y,z,t} \right)}} \\{= {f\left( {r,t} \right)}}\end{matrix}$where x,y,z are Cartesian coordinates in three dimensions; and

r is a three-dimensional positional vector.

The diffusive migration process can be expressed as a differentialequation, indicating the rate of change of concentration of malignantcells, c, with time at each point in a spatial domain over which thescalar concentration function c is defined. One approach is to model thediffusive migration process as:

$\frac{\partial{c\left( {r,t} \right)}}{\partial t} = {D{\nabla^{2}{c\left( {r,t} \right)}}}$

where

-   -   D is a constant diffusion coefficient over the entire domain of        the concentration function;    -   ∇² is the Laplacian operator; and    -   c(r, t) is a scalar function of the concentration of malignant        cells at position r at time t.        This provides a spherically symmetric diffusive migration,        which, as discussed above, generally does not occur in        anisotropic biological tissues. A second model is:

$\frac{\partial{c\left( {r,t} \right)}}{\partial t} = {\nabla{\cdot \left( {{D(r)}{\nabla{c\left( {r,t} \right)}}} \right)}}$

-   -   where D(r) is a scalar diffusion coefficient that varies over        the domain of the concentration function.        A third model employs a second-order-tensor field D to express        the rate of diffusion in three dimensions from each position r        within the domain of the concentration function:

$\frac{\partial c}{\partial t} = {\sum\limits_{i = 1}^{3}\;{\sum\limits_{j = 1}^{3}{\frac{\partial}{\partial x_{i}}\left( {{D_{ij}(r)}\frac{\partial c}{\partial x_{j}}} \right)}}}$where D(r) is a tensor field.

One computational model used in embodiments of the present invention fortumor growth is a differential equation with a cell-proliferationcomponent and a diffusive-migration component:

$\frac{\partial c}{\partial t} = {{\nabla{\cdot \left( {{D(r)}{\nabla c}} \right)}} + {\rho\;{c\left( {1 - \frac{c}{k}} \right)}}}$

where

-   -   D(r) is either a scalar function of position or a tensor field,        as discussed above;    -   k is the maximum concentration of malignant cells obtained in        the tissue; and    -   ρ is a rate of malignant-cell proliferation.        There is no closed-form solution to this equation. Instead,        numerical solutions are computationally generated by        numerical-analysis methods. This computational model is referred        to as the        diffusive-migration/bounded-exponential-cell-proliferation model        (“DM/BECP model”). The cell proliferation component is bounded        by the maximum-concentration constraint, discussed above, and        incorporated in the proliferation component. The rate of change        of the concentration of malignant cells with time may be        expressed in various different units, such as the number of        cells per mm³ per day, month, or year, with appropriate units        used for the diffusion coefficient, cell-proliferation constant,        and maximum malignant-cell concentration. The model is further        constrained by:        rεB        c(r,t=0)=c(r ₀)        n·∇c=0 on ∂B        These constraints include: (1) the position at which the        malignant-cell concentration is given by solution to the        differential equation is within the volume of a particular        tissue, B; (2) there is an initial malignant-cell concentration        c(r₀) at each point r; and (3) there is no migration of        malignant cells to points external from the volume of the        tissue.

The DM/BECP model is related to various other computational models usedto modeling various other problem domains. One such computational modelis the Fisher-Kolmogorov equation initially used to describe theadvancement of advantageous genes through a population. The DM/BECPmodel for tumor growth, the Fisher-Kolmogorov equation, and other suchequations are known to have a traveling-wave solution, in the case ofspherical spatial symmetry, in which a concentration wave expandsradially outward from an initial central point with a constant velocityν equal to:ν=2√{square root over (Dρ)}

FIGS. 6A-F illustrate tumor growth, in cross-sections through a tumorvolume, as monitored by two different types of medical imaging. Aglioblastoma is assumed in the example of FIGS. 6A-C. As discussedabove, the glioblastoma has a central, necrotic region 602 and anenclosing endemic region 604. The boundary of the internal, necroticregion is detected by T1-Gd-MRI (“T1”) imaging while the boundary of theendemic region 604 is detected by T2-MRI (“T2”) imaging. FIGS. 6A-C showthe growth of both regions within the tumor over time. Note that, ingeneral, diffusive invasion of tumor cells extends past theT2-imaging-detected boundary. FIG. 6D shows a plot of malignant-cellconcentration over time along a line x (606 in FIG. 6A) through thecross-section of the tumor shown in FIGS. 6A-C. The concentration curvefor FIG. 6A is shown as a solid line 610, the concentration curve forFIG. 6B is shown as a first dashed line 612, and the concentration curvefor FIG. 6C is shown as a second dashed line 614. The horizontal dashedline 616 shows the concentration threshold for T1-imaging detection andthe horizontal line 618 shows the concentration threshold for T2-imagingdetection, respectively. These concentrations are designated t₁ 620 andt₂ 622. As shown in FIG. 6D, the concentration curves at each point intime form a concentration-gradient wave front which moves outward inboth the positive and negative directions from the x=0 central point ofthe tumor 624. For any curve representing the concentration gradient ata particular time point, as shown in FIG. 6E, the x coordinate for thewave front, x_(WF), may be computed as the x coordinate of a point 630midway between the point of intersection of the concentration curve withthe horizontal t₁ and t₂ threshold concentration lines. Thepositive-wave-front coordinate, x_(WF), is plotted versus time, in FIG.6D, for tumor growth represented in FIGS. 6A-C as modeled by the DM/BECPmodel of the present invention. The position of the wave front varieslinearly with time, with the slope of the plotted line equal to theconstant wave-front velocity 2√{square root over (Dρ)}. In fact, thewave-front velocity may initially be non-linear, but approacheslinearity with increasing time.

The DM/BECP model, used in certain embodiments of the present invention,has been shown to accurately predict tumor growth for glioblastomas andother gliomas, and may, in addition, be applicable to a wide variety ofother types of tumors and growths within normal tissues. However, themodel is parameterized, with parameters ρ and D(r). These two parametersare tumor/patient specific. Thus, values for these parameters need to beobtained in order to use the DM/BECP model to characterize a particulartumor. As discussed above, assuming a constant diffusion coefficient D,the solution for the DM/BECP-model differential equation exhibits aconstant-velocity concentration wave front. Using two different medicalimages of a tumor, obtained at two different times t₁ and t₂, thevelocity of the concentration wave front is determined by:

$v = \frac{x_{{WF}{(t_{2})}} - x_{{WF}{(t_{1})}}}{t_{2} - t_{1}}$The wave fronts may be computed as equivalent radii of a spherical tumorof volume equal to the volumes of the tumor at times t₁ and t₂, using athreshold tumor cell concentration or other metric to determine thetumor boundary in the images. Alternatively, an average linear rate ofexpansion of the tumor along directions normal to surface of the tumoris computed as the wave-front velocity. Thus, two different MRI imagesare used to obtain one equation for the two unknown parameters D and ρ:

${2\sqrt{D\;\rho}} = \frac{x_{{WF}{(t_{2})}} - x_{{WF}{(t_{1})}}}{t_{2} - t_{1}}$

A second equation relating the two unknown parameters ρ and D is neededin order to determine values for the parameters ρ and D. An analyticalsolution for the differential equation

$\frac{\partial c}{\partial t} = {{D\;{\nabla^{2}c}} + {\rho\; c}}$in spherical symmetry with a point-source initial condition on aninfinite domain is known. From this analytical solution, an equation forthe ratio of ρ to D is derived as:

$\left( \frac{\rho}{D} \right)_{c} = {\left\lbrack \frac{4{\ln\left( \frac{t_{1}}{t_{2}} \right)}}{\left( {T\; 2r_{eqv}} \right)^{2} - \left( {T\; 1r_{eqv}} \right)^{2}} \right\rbrack\left\lbrack {\left\lbrack \frac{T\; 1r_{eqv}^{2}{\ln\left( \frac{t_{1}}{t_{2}} \right)}}{\left( {T\; 2r_{eqv}} \right)^{2} - \left( {T\; 1r_{eqv}} \right)^{2}} \right\rbrack + {\ln{\frac{t_{1}k}{N_{0}}\left\lbrack \frac{\Pi\left( {{T\; 2r_{eqv}^{2}} - {T\; 1r_{eqv}^{2}}} \right)}{\ln\left( \frac{t_{1}}{t_{2}} \right)} \right\rbrack}^{3/2}}} \right\rbrack}$

where

-   -   T1r_(eqv) and T2r_(eqv) are the radii of the equivalent,        spherical tumor volumes corresponding to the T1-imaging and        T2-imaging tumor volumes, respectively;    -   t₁ and t₂ are the threshold concentrations for T1-imaging and T2        imaging, respectively,    -   k is the maximum malignant-cell concentration; and    -   N₀ is an initial malignant-cell concentration or number.        The ratio

$\frac{\rho}{D}$computed from the above equation, with k and N₀ assumed to be constants,

${\left( \frac{\rho}{D} \right)_{c} = {f_{analytical}\left( {{T\; 1r_{eqv}},{T\; 2r_{eqv}}} \right)}},$is closely related to observed values for the ratio obtained fromsimulations using the DM/BECP model,

$\left( \frac{\rho}{D} \right)_{0} = {{model}\text{-}{simulation}\mspace{14mu}{\left( {{T\; 1r_{eqv}},{T\; 2r_{eqv}}} \right).}}$

FIGS. 7A-D illustrate preparation of a lookup table from whichcorrections to computed

$\frac{\rho}{D}$values are obtained according to one embodiment of the presentinvention. As shown in FIG. 7A, the analytical solution for thedifferential equation

$\frac{\partial c}{\partial t} = {{D\;{\nabla^{2}c}} + {\rho\; c}}$is used to derive an analytical function for computing

$\left( \frac{\rho}{D} \right)_{c},{\left( \frac{\rho}{D} \right)_{c} = {f_{analytical}\left( {{T\; 1r_{eqv}},{T\; 2r_{eqv}}} \right)}},$where T1r_(eqv) is the radius of a spherical volume equivalent to thevolume of the inner, necrotic portion of the glioblastoma, as measuredby T1-imaging, and T2r_(eqv) is the radius of the spherical volumeequivalent to the total tumor volume, as measured by T2-imaging. Theadditional parameters k and N₀ are treated as constants. Thus, as shownin FIG. 7A, a table of

$\left( \frac{\rho}{D} \right)_{c}$values are computed by the analytical function for various T1r_(eqv) andT2r_(eqv) values, and included in a table 702. Similarly, as shown inFIG. 7B, an additional table 704 of observed values

$\left( \frac{\rho}{D} \right)_{0}$are prepared by DM/BECP-model-based simulations of tumor growth, where

$\left( \frac{\rho}{D} \right)_{0} = {{model}\text{-}{simulation}\mspace{14mu}{\left( {{T\; 1r_{eqv}},{T\; 2r_{eqv}}} \right).}}$As shown in FIG. 7C, a table 706 of corrections, or Δ values,

${\Delta = {\left( \frac{\rho}{D} \right)_{0} - \left( \frac{\rho}{D} \right)_{c}}},$is prepared from tables 702 and 704 by each value in table 702 from thecorresponding value in table 704, as shown in FIG. 7C. Finally, as shownin FIG. 7D, a lookup table (“LUT”) is prepared from table 706, shown inFIG. 7C, by rearranging the table, so that, in the rearranged table, theΔ values are indexed by

$\left( \frac{\rho}{D} \right)_{c}$and T1r_(eqv).

Once the LUT (708 in FIG. 7D) is prepared, as discussed above, then anestimated observed

$\frac{\rho}{D}$value,

$\left( \frac{\rho}{D} \right)_{0},$is obtained from estimates of T1r_(eqv) and T2r_(eqv) computed from a T1and a T2 image, by:

$\left( \frac{\rho}{D} \right)_{0} = {{{LUT}\left( {\left( \frac{\rho}{D} \right)_{c},{T\; 1\; r_{eqv}}} \right)} + \left( \frac{\rho}{D} \right)_{c}}$This provides the second equation relating the two unknown parameters Dand ρ. With two equations for the two unknowns, values for ρ and D arecomputed by, for example:

$\frac{\rho}{D} = x_{1}$ ${2\sqrt{D\;\rho}} = x_{2}$${2\sqrt{D^{2}x_{1}}} = x_{2}$ $D = \frac{x_{2}}{2\sqrt{x_{1}}}$$\rho = \frac{x_{2}\sqrt{x_{1}}}{2}$

where x₁ and x₂ are numeric values computationally obtained from images,as discussed above.

With numerical values for D and ρ obtained as discussed above, acomplete computation model for a particular tumor/patient arenumerically generated. Thus, for example, a tumor volume can computedfor a malignant-cell concentration threshold t₃ lower than either t₁ ort₂ for performing various treatment plans. In addition, the ratio

$\frac{\rho}{D}$characterizes the aggressiveness of a particular tumor. Thus, the ratio

$\frac{\rho}{D}$can help to determine the prognosis for a patient, as well as to betterinform treatment selection and treatment planning. Other characteristicsof the tumor may be derived from various expressions derived from thediffusive-migration/exponential-cell-proliferation model.

FIG. 8 illustrates the complex, three-dimensional nature of the imagingand model-parameter-estimation problems associated with obtaining apatient-specific model for glioblastoma growth, according to embodimentsof the present invention. In FIG. 8, a Cartesian coordinate system 802is superimposed onto a cutaway illustration of a glioblastoma within apatent's cranium. In FIG. 8, three partial cross-sections 804-806 of aglioblastoma tumor are shown with respect to three mutually orthogonalplanes defined by the positive x, y, and z axes. The tumor has aconvoluted, irregular surface with many protrusions and invaginations.Those familiar with biological imaging and computing would immediatelyrecognize that even the relatively straightforward task of placing twodifferent images, taken at two different times, of the tumor into acommon coordinate system, given a generally irregular and complex tumorexpansion, represents an exceedingly tedious and difficult task that canfeasibly be undertaken only by computational methods. The spatialregistration task requires methodically exploring many differentrelative orientations of the two different images to determine arelative orientation that, in part, minimizes the spatial difference,computed in three dimensions, between two different tumor surfaces undervarious orientation constraints that may be available from imagingparameters. Such spatial registration tasks require millions ofcomputations. Other tasks, including computing tumor-expansionvelocities, ρ/D values, D tensor fields or D(x,y,z) scalar fields arealso computationally complex, and cannot be feasibly undertaken by anymeans other than computation on an electronic computer. It is notpossible to compute these characteristics and parameters accurately byhand, in any time frame, and quite impossible to carry out the millionsof calculations needed to compute these characteristics and parametersin reasonable and feasible time frames needed for medical prognosis,diagnosis, and treatment planning.

FIGS. 9A-E illustrate one method by which the scalar diffusion fieldD(r) is determined. The method illustrated in FIGS. 9A-E is but oneexample of various methods that are employed in order to determine theD(r) scalar field. As shown in FIG. 9A, two images of a tumor areobtained at different times t₁ and t₂. In FIG. 9A, the tumor is shown asdisk-shaped, idealized cross-sections, where an inner disk region 902represents a cross-section through the tumor volume at time t=t₁ anddisk-shaped region 904 represents a cross-section through the tumorvolume at time t=t₂. A three-dimensional grid, assumed in the currentdiscussion to be a lattice with cubic unit cells, is superimposed overthe region of tumor expansion 906 lying between the surface of the tumorat time t=t₁ 902 and the surface of the tumor at time t=t₂ 904, as shownin two-dimensional cross-section in FIG. 9B. Each three-dimensional cellwithin the three-dimensional grid is associated with two parameters κand γ, where:

$\kappa = \frac{{white}\text{-}{matter}\mspace{14mu}{cells}}{{total}\mspace{14mu}{cells}}$$\gamma = \frac{{gray}\text{-}{matter}\mspace{14mu}{cells}}{{total}\mspace{14mu}{cells}}$The relative numbers of white-matter cells κ and gray-matter cells γ foreach grid element is obtained from the t=t₁ image by variousimage-analysis techniques, or by use of the t=t₁ image along withgeneralized models of the distribution of white-matter and gray-mattercells within the brain.

In the determination of the scalar diffusion-coefficient field, variouspaths are constructed between points on the surface of the t=t₁ volumeand points on the surface of the t=t₂ volume. The paths are constructedas a sequence of adjacent grid elements, each pair of consecutive gridelements in a path sharing at least one common point in thethree-dimensional grid, with the first grid element in a path lying onthe t=t₁ surface and the final element in a grid-element path lying onthe t=t₂ surface. FIG. 9C shows, in cross-section, two differentpossible subpaths 910 and 912 leading from a first grid element 914 toone of two second grid elements 916 and 918, respectively. The length ofa subpath that includes grid elements 914 and 916 is equal to d_(g), thegrid spacing, whereas the length of a subpath that includes gridelements 914 and 918 is √{square root over (2)}d_(g). As shown in FIG.9D, there is a probability associated with each pair of grid elementsthat include a given grid element i 920 and each of the surrounding gridelements that share at least one point with grid element i for inclusionas a next subpath in a sequence of grid elements that represents a pathfrom the surface of the t=t₁ volume and the surface of the t=t₂ volume.In other words, at each point in the three-dimensional grid, aprobability P_(i,j) is calculated for a tumor surface to expand in anyof 26 different directions from grid element i to an adjacent gridelement j. Eight of the directions, shown in FIG. 9D, are planar, withan additional nine directions leading to grid elements below the planeof grid element i and an additional nine directions leading to gridelements above the plane of grid element i.

The method for determining the scalar diffusion-coefficient fieldassumes that the velocity of tumor expansion, at each point in time andsmall volume of space, is a linear combination of a velocity V_(w) inwhite matter and a velocity V_(g) in gray matter, with contributions ofvelocities V_(w) and V_(g) proportional to the relative amounts of whitematter and gray matter in the small volume:V=κV _(w) +γV _(g)As shown in FIG. 9E, the current method selects a large number of pointscoincident with the t=t₁ tumor surface, including points p₁₁, p₂₁ andp₃₁ 930-932 in FIG. 9E, and computes a most likely path to correspondingpoints on the t=t₂ surface of the tumor, shown in FIG. 9E as points p₁₂,p₂₂, and p₃₂ 934-936. The aggregate κ and γ values for a path p iscomputed as:

$\kappa_{p} = {\sum\limits_{i = 1}^{n}\kappa_{i}}$$\gamma_{p} = {\sum\limits_{i = 1}^{n}\gamma_{i}}$where there are n grid elements in the path. As discussed with referenceto FIG. 9C, path-component length is computed between each twosuccessive grid elements in a path as:

${{path\_ component}{\_ length}} = \left\{ \begin{matrix}{{{when}\mspace{14mu}{planar}\mspace{14mu}{diagonal}},} & \sqrt{2\; d_{g}} \\{{{when}\mspace{14mu}{non}\text{-}{planar}\mspace{14mu}{diagnonal}},} & {1.73\; d_{g}} \\{otherwise} & d_{g}\end{matrix} \right.$The length of a path p is then:

${{length}(p)} = {\sum\limits_{i = 1}^{n - 1}{{path\_ component}{{\_ length}_{i}.}}}$The velocity of tumor expansion along a given path p, V_(p), is computedas:

$V_{p} = \frac{{length}(p)}{t_{2} - t_{1}}$From any grid element j₁ lying on the t=t₁ surface, there is acomputable probability P(j₁, j₂, . . . j_(n)) that the tumor expandsfrom grid element j₁ to grid element j_(n) on the t=t₂ surface along thesequence of grid elements j₁, j₂, . . . j_(n):P(j ₁ ,j ₂ , . . . ,j _(n))=f(P _(j) ₁ _(,j) ₂ ,P _(j) ₂ _(,j) ₃ , . . .,P _(j) _(n−1) _(,j) _(n) )The function ƒ( ) is computed in various ways, including combining theprobabilities of tumor expansion between adjacent grid-element pathcomponents in ways that consider the κ and γ values for the gridelements in the path, the distance of the grid elements from t=t₁surface, and other parameters. Assuming a computable function ƒ( ), thepaths between each grid element m on the t=t₁ surface to a correspondinggrid element n on the t=t₂ surface is computed as the path with maximumcomputed probability:

$P_{m,n} = {\max\limits_{{{{all}\mspace{14mu}{possible}\mspace{14mu} j_{1}} = m},j_{2},\ldots\mspace{14mu},{j_{n} = n}}{P\left( {}_{j_{1},j_{2},\ldots\mspace{14mu},j_{n}} \right)}}$After computing paths from some well-distributed set of t=t₁ surfacepoints to the t=t₂ surface, a set of paths p₁, p₂, . . . p_(n) isobtained. The velocities of tumor expansion for each path are thencomputed as:

${\begin{bmatrix}\kappa_{1} & \gamma_{1} \\\kappa_{2} & \gamma_{2} \\\vdots & \vdots \\\kappa_{n - 1} & \gamma_{n - 1} \\\kappa_{n} & \gamma_{n}\end{bmatrix}\begin{bmatrix}V_{w} \\V_{g}\end{bmatrix}} = \begin{bmatrix}V_{1} \\\vdots \\V_{n - 1} \\V_{n}\end{bmatrix}$This is an overdetermined set of equations that allows the velocitiesV_(w) and V_(g) to be computed by a least-squares regression or by othertechniques for determining variable values in overdetermined systems ofequations. Finally, the diffusion coefficients D_(w) and D_(g) areobtained by assuming a constant proliferation coefficient ρ as follows:

$D_{w} = \frac{V_{w}^{2}}{4\rho}$ $D_{g} = \frac{V_{g}^{2}}{4\rho}$For each grid element, a scalar diffusion coefficient is computed as anumeric combination of the diffusion coefficients D_(w) and D_(g) basedon the κ and γ values for the grid element.

The DM/BECP model can be further expanded to take into account not onlythe heterogeneity of the brain tissue, but also its anisotropy,including facilitation of glioma cell migration in the direction ofwhite matter fibers, as briefly mentioned above. The model is expressedas:

${\frac{\partial c}{\partial t} = {{\nabla{\cdot \left( {{D(x)}{\nabla c}} \right)}} + {f(c)}}},$where D is a diffusion tensor, at each point within a tissue, thatdescribes tumor cell diffusion. The diffusion tensor is a 3-by-3symmetric positive definite matrix that models the local anisotropy ofthe cell diffusion tensor. The expanded model takes into account bothlocation and direction within the structure of the brain tissue, withproliferation ƒ(c) modeled by a logistic law with c_(m)=10⁵ cells/mm³.

Initial conditions are defined as c(x, t=0)=c₀(x) and, to inhibitmigration of glioma cells outside of the brain tissue, D(x)∇c·n=0 for xon the sulcal and ventricular boundary of the brain, where n is thesurface normal. The above model equation can be solved, to estimatetumor-cell density c, using 3D finite differences, which involvesapproximating temporal and spatial derivatives by their discreteexpressions. First, the following PDE is solved:

$\quad\left\{ \begin{matrix}{\frac{\partial c}{\partial t} = {{\nabla{\cdot \left( {{D(x)}{\nabla c}} \right)}} + {\rho\; c}}} \\{{c\left( {x,{t = 0}} \right)} = {c_{0}(x)}}\end{matrix} \right.$where c₀ is the concentration of the tumor cells at t=0. Each part ofthis PDE is then discretized using ΔT as a time step and (ΔX,ΔY,ΔZ) asspace steps to obtain the following discrete equation,

${\frac{C^{n + 1} - C^{n}}{\Delta\; T} = {AC}^{n}},$where C^(n) is a vector representing cell density at the time nΔT and Ais a large sparse matrix representing discrete operators. This discreteequation represents the forward Euler method. Instead of using theforward Euler method, an alternative θ-method is used in order toincrease numerical stability. The θ-method consists of combining forwardand backward numerical schemes, as follows:

${\frac{C^{n + 1} - C^{n}}{\Delta\; T} = {{\left( {1 - \theta} \right){AC}^{n}} + {\theta\;{AC}^{n + 1}}}},$

Denoting δ^(n)=C^(n+1)−C^(n), the cell concentration at time (n+1)ΔT iscomputed by first solving the equation(I−θΔTA)δ^(n) =ΔTAC ^(n),where I is the identity matrix, which gives the value for δ^(n). Thepreconstrained gradient method is used to solve this large sparse systemand obtain the tumor-cell concentration asC ^(n+1) =C ^(n)+δ^(n)Initial conditions are represented by a tumor-cell concentration C⁰ in amanually-selected grid element. Using the notation(λ_(i),e_(i))_(i=1,2,3) to denote the eigenvalues and eigenvectors ofthe water diffusion tensor D obtained by diffusion-tensor imaging(“DTI”). The following transformed tensor D is to represent thediffusion tensor of tumor cells:D= λ ₁(r)e ₁ e ₁ ^(T)+ λ ₂(r)e ₂ e ₂ ^(T)+ λ ₃(r)e ₃ e ₃ ^(T).This operation modifies the eigenvalues of the tensor, but not theeigenvectors, which implies that the tensor orientation is preservedwhile the diffusion values along the principal axes and the anisotropyare changed. The transformation is controlled by a parameter r. When requals 1, the tensor is not changed. r<1 corresponds to a decrease inanisotropy and r>1 to an increase in anisotropy.

The finite-differences method is next described. Let c(t,x) denote thecell concentration at time t and location x of the space. The expandedmodel equation can be expressed as:

${\frac{\partial c}{\partial t} = {{\sum\limits_{i,{j = 1}}^{3}\;{D_{ij}\frac{\partial^{2}c}{{\partial x_{i}}{\partial x_{j}}}}} + {\sum\limits_{i = 1}^{3}{{\overset{\sim}{D}}_{i}\frac{\partial c}{\partial x_{i}}}} + {\rho\; c}}},{where}$${\overset{\sim}{D}}_{i} = {{\left( {D\nabla} \right)i} = {\sum\limits_{j = 1}^{3}{\frac{\partial D_{ij}}{\partial x_{j}}.}}}${tilde over (D)} is derived from the tensor D by computing itselement-by-element gradient. Let c(nΔT,iΔX,jΔY,kΔZ)≡C_(i,j,k) ^(n) andlet C^(n) denote the vectorized version of C_(i,j,k) ^(n). The followingdiscrete schemes are used for computing the derivatives:

$\left. \frac{\partial c}{\partial t}\rightarrow\frac{C_{i,j,k}^{n + 1} - C_{i,j,k}^{n}}{\Delta\; T} \right.$$\left. \frac{\partial c}{\partial x}\rightarrow\frac{C_{{i + 1},j,k}^{n} - C_{{i - 1},j,k}^{n}}{2\Delta\; X} \right.$$\left. \frac{\partial^{2}c}{\partial x^{2}}\rightarrow\frac{C_{{i + 1},j,k}^{n} - {2C_{i,j,k}^{n}} + C_{{i - 1},i,j,k}^{n}}{\Delta\; X^{2}} \right.$$\left. \frac{\partial^{2}c}{{\partial x}{\partial y}}\rightarrow\frac{C_{{i + 1},{j + 1},k}^{n} + C_{{i - 1},{j - 1},k}^{n} - C_{{i + 1},{j - 1},k}^{n} - C_{{i - 1},{j + 1},k}^{n}}{4\Delta\; X\;\Delta\; Y} \right.$and similarly for

$\frac{\partial c}{\partial y},\frac{\partial c}{\partial z},\frac{\partial^{2}c}{\partial y^{2}},\frac{\partial^{2}c}{\partial z^{2}},\frac{\partial^{2}c}{{\partial x}{\partial z}},\;{{and}\mspace{14mu}{\frac{\partial^{2}c}{{\partial y}{\partial z}}.}}$

The following discrete model-equation expression for is then obtained:

$\frac{C_{i,j,k}^{n + 1} - C_{i,j,k}^{n}}{\Delta\; T} = {{D_{11}\frac{C_{{i + 1},j,k}^{n} - {2C_{i,j,k}^{n}} + C_{{i - 1},j,k}^{n}}{\Delta\; X^{2}}} + {D_{22}\frac{C_{i,{j + 1},k}^{n} - {2C_{i,j,k}^{n}} + C_{i,{j - 1},k}^{n}}{\Delta\; Y^{2}}} + {D_{33}\frac{C_{i,j,{k + 1}}^{n} - {2C_{i,j,k}^{n}} + C_{i,j,{k - 1}}^{n}}{\Delta\; Z^{2}}} + {2\; D_{12}\frac{C_{{i + 1},{j + 1},k}^{n} + C_{{i - 1},{j - 1},k}^{n} - \left( {C_{{i + 1},{j - 1},k}^{n} + C_{{i - 1},{j + 1},k}^{n}} \right)}{4\Delta\; X\;\Delta\; Y}} + {2\; D_{13}\frac{C_{{i + 1},j,{k + 1}}^{n} + C_{{i - 1},j,{k - 1}}^{n} - \left( {C_{{i + 1},j,{k - 1}}^{n} + C_{{i - 1},j,{k + 1}}} \right)}{4\Delta\; X\;\Delta\; Z}} + {2\; D_{23}\frac{C_{i,{j + 1},{k + 1}}^{n} + C_{i,{j - 1},{k - 1}}^{n} - \left( {C_{i,{j + 1},{k - 1}}^{n} + C_{i,{j - 1},{k + 1}}} \right)}{4\Delta\; Y\;\Delta\; Z}} + {{\overset{\sim}{D}}_{1}\frac{C_{{i + 1},j,k}^{n} - C_{{i - 1},j,k}^{n}}{2\Delta\; X}} + {{\overset{\sim}{D}}_{2}\frac{C_{i,{j + 1},k}^{n} - C_{i,{j - 1},k}^{n}}{2\Delta\; Y}} + {{\overset{\sim}{D}}_{3}\frac{C_{i,j,{k + 1}}^{n} - C_{i,j,{k - 1}}^{n}}{2\Delta\; Z}} + {\rho\; C_{i,j,k}^{n}}}$which can be rewritten as

${\frac{C^{n + 1} - C^{n}}{\Delta\; T} = {AC}^{n}},$where A is a big sparse matrix with 19 diagonal elements and whereelements of A depend on D, {tilde over (D)}, ρ, and on the space/timediscretization step sizes (ΔT,ΔX,ΔY,ΔZ).

To fit the best solution, the parameters are chosen ΔT and θ with care.Neglecting the diffusion term, the following recursive scheme isobtained:

$\frac{C^{n + 1} - C^{n}}{\Delta\; T} = {{\left( {1 - \theta} \right)\rho\; C^{n}} + {\theta\;\rho\; C^{n + 1}}}$$C^{n + 1} = {\frac{1 - {{\rho\Delta}\; T\;\theta} + {{\rho\Delta}\; T}}{1 - {{\rho\Delta}\; T\;\theta}}C^{n}}$Denoting u=ρΔT, the following solution for an initial condition C⁰ isobtained:

$C^{n} = {\left( \frac{1 - {\theta\; u} + u}{1 - {\theta\; u}} \right)^{n}{C^{0}.}}$The analytical solution for this exponential proliferation law isC^(n)=C⁰e^(nρΔT), leading to the relation

${\left( \frac{1 - {\theta\; u} + u}{1 - {\theta\; u}} \right)^{n} = {\mathbb{e}}^{nu}},$which gives a relationship between the parameters

$\theta = {{\frac{1}{u} - \frac{1}{{\mathbb{e}}^{u} - 1}} = {\frac{1}{{\rho\Delta}\; T} - {\frac{1}{{\mathbb{e}}^{{\rho\Delta}\; T} - 1}.}}}$In this case, the diffusion D is zero.

The DM/BECP model, discussed above, can be expanded in variousadditional ways. As one example, the model can be expanded to include anadditional term for the effects of radiation therapy:

$\frac{\partial c}{\partial t} = {{\nabla{\cdot \left( {D{\nabla c}} \right)}} + {\rho\;{c\left( {1 - \frac{c}{k}} \right)}} - {{Rc}\left( {1 - \frac{c}{k}} \right)}}$where

${R\left( {x,t,{{Dose}\left( {x,t} \right)}} \right)} \equiv \left\{ \begin{matrix}0 & {{{for}\mspace{14mu} t} \notin {therapy}} \\{1 - {S\left( {\alpha,\beta,{{Dose}\left( {x,t} \right)}} \right)}} & {{{for}\mspace{14mu} t} \in {therapy}}\end{matrix} \right.$In this equation, the term S(α, β, Dose(x,t)) represents the survivalrate of cells exposed to a particular dose of radiation, where α and βare parameters for a well-known radiation-exposure model.

As another example, the model can be expanded to take into accountangiogenesis within a tumor. As summarized below, the expandedcomputational model that includes the angiogenic cascade, referred to asthe “PIHNA model,” is expressed as five linked partial differentialequations describing the interactions between normoxic (c), hypoxic (h),necrotic (n) and endothelial cells (ν) with angiogenic factors (a):

$\frac{\partial c}{\partial t} = {\overset{\overset{\begin{matrix}{{Net}\mspace{14mu}{dispersal}\mspace{14mu}{of}} \\{normoxic} \\{{glioma}\mspace{14mu}{cells}}\end{matrix}}{︷}}{\nabla{\cdot \left( {{D\left( {1 - T} \right)}{\nabla c}} \right)}} + \overset{\overset{\begin{matrix}{{Net}\mspace{14mu}{proliferation}\mspace{14mu}{of}} \\{normoxic} \\{{glioma}\mspace{14mu}{cells}}\end{matrix}}{︷}}{\rho\;{c\left( {1 - T} \right)}} + \overset{\overset{\begin{matrix}{{Conversion}\mspace{14mu}{of}} \\{{hypoxic}\mspace{14mu}{to}} \\{normoxic}\end{matrix}}{︷}}{\gamma\;{hV}} - \overset{\overset{\begin{matrix}{Conversion} \\{{of}\mspace{14mu}{normoxic}\mspace{14mu}{to}} \\{hypoxic}\end{matrix}}{︷}}{\beta\;{c\left( {1 - V} \right)}} - \overset{\overset{\begin{matrix}{Conversion} \\{{of}\mspace{14mu}{normoxic}\mspace{14mu}{to}} \\{necrotic}\end{matrix}}{︷}}{\alpha_{n}{nc}}}$$\frac{\partial h}{\partial t} = {\overset{\overset{\begin{matrix}{{{Dispersa}l}\mspace{14mu}{of}} \\{hypoxic} \\{{glioma}\mspace{14mu}{cells}}\end{matrix}}{︷}}{\nabla{\cdot \left( {{D\left( {1 - T} \right)}{\nabla h}} \right)}} - \overset{\overset{\begin{matrix}{{Conversion}\mspace{14mu}{of}} \\{{hypoxic}\mspace{14mu}{to}} \\{normoxic}\end{matrix}}{︷}}{\gamma\;{hV}} + \overset{\overset{\begin{matrix}{Conversion} \\{{of}\mspace{14mu}{normoxic}\mspace{14mu}{to}} \\{hypoxic}\end{matrix}}{︷}}{\beta\;{c\left( {1 - V} \right)}} - \overset{\overset{\begin{matrix}{Conversion} \\{{of}\mspace{14mu}{hypoxicc}\mspace{20mu}{to}} \\{necrotic}\end{matrix}}{︷}}{\left( {{\alpha_{h}{h\left( {1 - V} \right)}} + {\alpha_{n}{nh}}} \right)}}$$\mspace{79mu}{\frac{\partial n}{\partial t} = \overset{\overset{{{Conversion}\mspace{14mu}{of}\mspace{14mu}{hypoxic}},\mspace{11mu}{{normoxic}\mspace{14mu}{and}\mspace{14mu}{vasculature}\mspace{20mu}{to}\mspace{14mu}{necrotic}}}{︷}}{{\alpha_{h}{h\left( {1 - V} \right)}} + {\alpha_{n}{n\left( {c + h + v} \right)}}}}$$\mspace{79mu}{\frac{\partial v}{\partial t} = {\overset{\overset{\begin{matrix}{{Dispersal}\mspace{14mu}{of}} \\{vasculature}\end{matrix}}{︷}}{\nabla{\cdot \left( {{D_{v}\left( {1 - T} \right)}{\nabla v}} \right)}} + \overset{\overset{\begin{matrix}{{Net}\mspace{14mu}{proliferation}} \\{{of}\mspace{14mu}{vasculature}}\end{matrix}}{︷}}{\mu\frac{a}{K_{m} + a}{v\left( {1 - T} \right)}} - \overset{\overset{\begin{matrix}{{Conversion}\mspace{14mu}{of}\mspace{14mu}{vasculature}} \\{{to}\mspace{14mu}{necrotic}}\end{matrix}}{︷}}{\alpha_{n}{nv}}}}$$\frac{\partial a}{\partial t} = {\overset{\overset{\begin{matrix}{{Diffusion}\mspace{14mu}{of}} \\{angiogenic} \\{factors}\end{matrix}}{︷}}{\nabla{\cdot \left( {D_{a}{\nabla a}} \right)}} + \overset{\overset{\begin{matrix}{{Net}\mspace{14mu}{production}} \\{{of}\mspace{14mu}{angiogenic}\mspace{14mu}{factors}}\end{matrix}}{︷}}{{\delta_{c}c} + {\delta_{h}h}} - \overset{\overset{\begin{matrix}{{Net}\mspace{14mu}{consumption}\mspace{14mu}{of}} \\{{angiogenic}\mspace{14mu}{factors}}\end{matrix}}{︷}}{q\;\mu\frac{a}{K_{m} + a}{v\left( {1 - T} \right)}} - {\omega\;{av}} - \overset{\overset{\begin{matrix}{{Decay}\mspace{14mu}{of}} \\{{angiogenic}\mspace{14mu}{factors}}\end{matrix}}{︷}}{\lambda\; a}}$where V=ν/(c+h+ν) and T=(c+h+ν+n)/K.

The PIHNA-model definition begins with normoxic (c) glioma cells thatare assumed to move randomly away from high tumor cell concentration(i.e., diffuse) in the tissue at a rate D and to proliferate at a rateρ. These normoxic glioma cells compete for space with other neighboringcells. When the oxygen supply, determined by the relative amount ofvasculature supplying the tissue, is insufficient, normoxic glioma cellsconvert to a hypoxic state h at a rate β. When oxygen continues to beinadequate, the hypoxic cells become necrotic at a rate α. Both of theseprocesses correlate positively with the metabolic demands of the tumor,so β and α_(h) are proportional to ρ. On the other hand, when the oxygensupply increases at a later time and becomes sufficient, then hypoxicglioma cells will convert back to a normoxic state at a rate γ. It ispossible for γ to correlate negatively with ρ, again by reasoning ofmetabolism; however, γ is assumed fixed, with α_(h) considered toadequately describe the greater tendency towards hypoxia caused byincreased metabolic demands. Necrotic cells that come into contact withother cells can induce those viable cells to also become necrotic at arate α_(n).

FIG. 10 illustrates an alternative method for estimating the parametersρ and D according to one embodiment of the present invention. As shownin FIG. 10, a matrix of tumor-growth models, with specific values ofparameters ρ and D, is constructed using the above-discussed PIHNAmodel. In FIG. 10, this matrix of models is represented as atwo-dimensional plane 1002 indexed by a ρ-value axis 1004 and a D-valueaxis 1006. For each value of

$\frac{\rho}{D},$there is a line, such as line 1006, with slope

$\frac{\rho}{D}$extending from the origin 1010 outward. In the first-described two-imagemethod, the velocity of the tumor-expansion wavefront is computed, andthe velocity dependence on ρ and D is used as a second equation in twounknowns to allow ρ and D to be computed. However, the equivalent radiusof the necrotic tissue within a tumor, T0, is obtained from T1 and/or T2MRI images obtained at a single point in time, to provide a thirdobserved value, or parameter. It turns out that the equivalent radius ofthe tumor, T0, is computed from the PIHNA models created for eachdiscrete ρ and D value, represented by the two-dimensional plane 1002 inFIG. 10. The models that predict particular T0 values form isocurves,such as isocurve 1012, in the ρ-and-D-parameterized model space. Theratio

$\frac{\rho}{D}$can be computed from a T1 and a T2 image as discussed above withreference to the first-described two-image method, by an analyticalmethod. The equivalent radius of the necrotic portion of the tumor T0can be computed from either the T1 or T2 image, or by one or moreadditional types of imaging techniques. In order to determine ρ and D,the intersection between the line with slope

$\frac{\rho}{D},$such as line 1006 in FIG. 10, and the isocurve corresponding to theobserved equivalent radius of the necrotic region of the tumor, such asisocurve 1014 in FIG. 10, is determined. In FIG. 10, an intersectionbetween

$\frac{\rho}{D}$line 1006 and isocurve 1014 is shown as point 1016. The ρ and D valuescan then be directly obtained from the intersection point as the ρ and Dcoordinates of the intersection point, illustrated by dashed lines 1020and 1022 and the intersection points 1024 and 1026 of the dashed linesand the coordinate axes in FIG. 10. Therefore, ρ and D parameters thatcharacterize a tumor is obtained, by the first-discussed method, fromtwo different MRI images taken at two different times, or, using theabove-described alternative method, from a single T1, T2 image and froma pre-computed ρ-and-D-parameterized space of PIHNA computational modelsor other computational models for tumor growth. In practice, rather thanobtaining a single intersection point, a small region of intersection isobtained, reflecting imprecision in estimating T0 and

$\frac{\rho}{D},$providing generally small ranges of ρ and D values.

Extensive analysis and characterization of glioblastomas in glioblastomapatients has revealed that the parameters ρ and D, corresponding toproliferation and invasion parameters for tumor growth, do, indeed,characterize clinically observed glioblastoma growth in glioblastomapatients. The ratio

$\frac{\rho}{D}$can be interpreted as the relative steepness of the gradient of invadingcells peripheral to the tumor surface as seen in a T1 image. Thus, theratio

$\frac{\rho}{D}$is proportional to the proliferative nature of a tumor. Tumors with high

$\frac{\rho}{D}$ratios would be expected to show distinct borders and elevated necrosisassociated with hypoxia due to crowding and other microenvironmentalconstraints inherent to highly proliferative tumors with relatively lowinvasion. The ratio

$\frac{D}{\rho},$often referred to as the “invisibility index,” provides an estimate ofthe number or concentration of diffusely invading tumor cells peripheralto the imaged tumor boundary. Tumors with high

$\frac{D}{\rho}$values generally require applying more aggressive treatment to acontaining volume significantly larger than the observed tumor volume.The determination of ρ and D parameter values, using biological imagingcombined with computational modeling, thus provides useful andeasy-to-interpret quantitative characteristics of tumor growth that candirectly inform medical prognosis and treatment planning.

FIG. 11 provides a control-flow diagram for a routine “tumor-growthanalysis” that represents one embodiment of the present invention. Thismethod is carried out in an electronic computing system under softwarecontrol. In step 1102, an initial medical image of the tumor is obtainedat time t₁. In step 1104, a second image of the tumor is obtained attime t₂, using the same imaging technique. In step 1106, either acomplementary image, using a different imaging technique, is obtained attime t₂ or a complementary image pair is obtained at time t₃. In step1108, a first numerical value is calculated for the wave-front velocityν, as discussed above, using the initial and second images. In step1110, a second numerical value for the ratio

$\frac{\rho}{D}$is computed using an analytical function and a lookup table, asdiscussed above. In step 1112, D and ρ are computed from the first andsecond numerical values, as discussed above. In addition, a scalar ortensor field D(r) for a spatially varying diffusion coefficient can bedetermined, using various methods discussed above. A numericalDM/BECP-model is then obtained, in step 1114. In step 1116, theparameters D, ρ, and, optionally, D(r) may be stored in an electronicdatabase for the patient, so that the DM/BECP-model can be subsequentlyregenerated. In step 1118, the routine “use model” is called to providevarious derived characteristics for the tumor.

FIG. 12 provides a control-flow diagram for a routine “use model” thatrepresents one embodiment of the present invention. In step 1202, theroutine “use model” determines whether a particular DM/BECP model iscurrently in memory, and, when not, uses a patient identifier toretrieve the DM/BECP model from an electronic database. The entireDM/BECP model may be retrieved or, alternatively, the DM/BECP model maybe regenerated from retrieved parameter values. Then, in the for-loop ofsteps 1207-1214, requests for results derived from the DM/BECP model orfor parameter values of the DM/BECP model are received from a user, andthe requested derived results and parameter values output to a user orto an electronic data-processing or data-storage device on behalf of auser. Derived results may include a representation of the spatial extentof the tumor, in three dimensions, based on a threshold malignant-cellconcentration, the

$\frac{\rho}{D}$ratio, indicative of tumor aggressiveness, and any of many other derivedresults. These results may be then used, by medical personnel, fordisease prognosis, for therapy planning, for data-collection purposes tofurther medical research, and for a variety of additional purposes.

FIG. 13 provides a control-flow diagram for the alternative,single-image method for determining ρ and D parameters for a tumor. Instep 1302, T1 and T2 MRI images are obtained from a patient at a singlepoint in time. In step 1304, the ratio

$\frac{\rho}{D}$is computed by the analytical method discussed above with reference tothe two-image method. In step 1306, the equivalent radius of thenecrotic region of the tumor, T0, is computed from the T1, T2 MRI image.In step 1308, the intersection point between a line with slope

$\frac{\rho}{D}$and the isocurve for the computed T0 value is found inρ-and-D-parameterized model space, allowing ρ and D to be obtained, instep 1310, as the ρ-and-D coordinates of the intersection point betweenthe line and isocurve.

Although the present invention has been described in terms of particularembodiments, it is not intended that the invention be limited to theseembodiments. Modifications will be apparent to those skilled in the art.For example, models other than the DM/BECP model may be employed tomodel various types of tumors, the models containing parametersdifferent from, or in addition to, D and ρ. Imaging techniques, otherthan T1 imaging and T2 imaging, may be used to determine tumor volumesand regions of hypointensity, of T0 regions, used for computing D and ρ.Additional types of imaging, using labeled DNA monomer analogues,labeled metabolites, and other substances, can be used to obtainadditional characterization of tumor volumes, necrotic volumes withintumors, and additional parameters and characteristics of tumors that canbe incorporated into additionally expanded computational models. In theabove discussion, a first, two-time-point based method for obtaining Dand ρ was discussed. This method is easily extended to makes use ofadditional data generated at additional time points, and can thus beconsidered to be an n-time-point method. Additional data can providemore accurate estimation of D and ρ in the general, n-time-point method.A variety of implementations of the present invention may be obtained byvarying development parameters, including programming language,operating system platform, control structures, data structures, modularorganization, and other such parameters. Method embodiments of thepresent invention may be incorporated into imaging and/or image-analysisequipment, as well as in various hospital information systems.

The foregoing description, for purposes of explanation, used specificnomenclature to provide a thorough understanding of the invention.However, it will be apparent to one skilled in the art that the specificdetails are not required in order to practice the invention. Theforegoing descriptions of specific embodiments of the present inventionare presented for purpose of illustration and description. They are notintended to be exhaustive or to limit the invention to the precise formsdisclosed. Many modifications and variations are possible in view of theabove teachings. The embodiments are shown and described in order tobest explain the principles of the invention and its practicalapplications, to thereby enable others skilled in the art to bestutilize the invention and various embodiments with various modificationsas are suited to the particular use contemplated.

It is intended that the scope of the invention be defined by thefollowing claims and their equivalents:
 1. A computer-implemented methodof modeling tumor growth for guiding medical intervention, comprising:(a) calculating ρ/D, an observed ratio of a malignant cell proliferationcoefficient, ρ, to a diffusion coefficient, D, from one or moredifferent measures of tumor volume from either a single or multipleimages of a tumor, which calculating comprises mapping the followingformula derived from an analytic solution of a linear tumor model to anon-linear tumor model:$\left( \frac{\rho}{D} \right)_{c} = {\left\lbrack \frac{4\;{\ln\left( \frac{t_{1}}{t_{2}} \right)}}{\left( {T\; 2r_{eqv}} \right)^{2} - \left( {T\; 1r_{eqv}} \right)^{2}} \right\rbrack\left\lbrack {\left\lbrack \frac{T\; 1r_{eqv}^{2}{\ln\left( \frac{t_{1}}{t_{2}} \right)}}{\left( {T\; 2r_{eqv}} \right)^{2} - \left( {T\; 1\; r_{eqv}} \right)^{2}} \right\rbrack + {\ln\;{\frac{t_{1}k}{N_{0}}\left\lbrack \frac{\prod\;\left( {{T\; 2r_{eqv}^{2}} - {T\; 1r_{eqv}^{2}}} \right)}{\ln\left( \frac{t_{1}}{t_{2}} \right)} \right\rbrack}^{3/2}}} \right\rbrack}$where T1r_(eqv) and T2r_(eqv) are the radii of the equivalent, sphericaltumor volumes corresponding to the one or more different measures oftumor volume, respectively; t₁ and t₂ are threshold tumor cellconcentrations for two of said measures of tumor volume, k is maximummalignant-cell concentration; and N₀ is an initial malignant-cellconcentration or number; (b) obtaining a value for D and a value for ρ;and (c) generating a computational model for growth of the tumor in acancer patient using the values of D and ρ obtained in step (b); whereinsteps (a)-(c) are performed using a computer system.
 2. The method ofclaim 1, further comprising calculating a tensor diffusion coefficient,D(r), prior to step (c), to generate the computational model to generatethe computational model.
 3. The method of claim 2, further comprisingstoring the computational model, D, ρ, D(r) or a combination thereof inan electronic database.
 4. The method of claim 1, wherein the measuresof tumor volume are from two different types of medical imagingperformed at approximately the same time.
 5. The method of claim 1,wherein, given a size of the tumor necrotic region (T0) measured from animage, the values of D and ρ are obtained by finding an intersectionpoint between: (i) a line with slope ρ/D; and (ii) an isocurvecorresponding to the observed equivalent radius of T0 provided by anon-linear proliferation invasion hypoxia necrosis angiogenesis (PIHNA)model isocurve.
 6. The method of claim 1, wherein a first and a secondimage are obtained (i) at two separate times by the same imagingtechnique; or (ii) at approximately the same time using differentimaging techniques.
 7. The method of claim 1, wherein aggressivenessand/or prognosis of the tumor is determined by ρ/D.
 8. The method ofclaim 1, wherein the computational model comprises one or moreadditional terms to predict the effects of cancer therapy, to modelangiogenesis within the tumor, or both.
 9. The method of claim 1,wherein the tumor is a glioma.
 10. The method of claim 1, furthercomprising calculating a concentration wave-front velocity of tumorexpansion using one or more of said images.
 11. The method of claim 10,wherein ρ/D is calculated from an image or images obtained at twoseparate times, time₀ and time₁, using the same imaging technique. 12.The method of claim 11, wherein the concentration wavefront velocity oftumor expansion is calculated as (i) equivalent radii of a sphericaltumor of volume equal to the volumes of the tumor at times time₀ andtime₁ using a threshold tumor cell concentration; or as (ii) an averagelinear rate of expansion of the tumor along directions normal to a thetumor's surface.
 13. The method of claim 10, wherein the observed ratioρ/D is obtained from a lookup table.
 14. The method of claim 13, whereinthe values for D and ρ are obtained using the concentration wave-frontvelocity and observed ratio of ρ/D.
 15. A computer program product formodeling tumor growth and for guiding medical intervention comprising anon-transitory computer readable medium which, when executed by acomputer causes a method to be performed, the method comprising: (a)calculating a ratio, ρ/D, an observed ratio of a malignant cellproliferation coefficient, ρ, to a diffusion coefficient, D, from one ormore different measures of tumor volume from either a single or multipleimages of a tumor, which calculating comprises mapping the followingformula derived from an analytic solution of a linear tumor model to anon-linear tumor model:$\left( \frac{\rho}{D} \right)_{c} = {\left\lbrack \frac{4\;{\ln\left( \frac{t_{1}}{t_{2}} \right)}}{\left( {T\; 2r_{eqv}} \right)^{2} - \left( {T\; 1r_{eqv}} \right)^{2}} \right\rbrack\left\lbrack {\left\lbrack \frac{T\; 1r_{eqv}^{2}{\ln\left( \frac{t_{1}}{t_{2}} \right)}}{\left( {T\; 2r_{eqv}} \right)^{2} - \left( {T\; 1\; r_{eqv}} \right)^{2}} \right\rbrack + {\ln\;{\frac{t_{1}k}{N_{0}}\left\lbrack \frac{\prod\;\left( {{T\; 2r_{eqv}^{2}} - {T\; 1r_{eqv}^{2}}} \right)}{\ln\left( \frac{t_{1}}{t_{2}} \right)} \right\rbrack}^{3/2}}} \right\rbrack}$where T1r_(eqv) and T2r_(eqv) are the radii of the equivalent, sphericaltumor volumes corresponding to the one or more different measures oftumor volume, respectively; t₁ and t₂ are threshold tumor cellconcentrations for two of said measures of tumor volume, k is maximummalignant-cell concentration; and N₀ is an initial malignant-cellconcentration or number; (b) obtaining a value for D and a value for ρ;and (c) generating a computational model for tumor growth of the tumorin a cancer patient using the values of D and ρ obtained in (b); whereinthe computer-readable medium is non-transitory.
 16. The computer programproduct of claim 15, wherein the method further comprises calculating atensor diffusion coefficient, D(r), prior to step (c) to generate thecomputational model.
 17. The computer program product of claim 16,wherein the method further comprises storing the computational model, D,ρ, D(r) or a combination thereof in an electronic database.
 18. Thecomputer program product of claim 16, wherein the measures of tumorvolume are performed by two different types of medical imaging taken atapproximately the same time.
 19. The computer program product of claim15, wherein, given a size of the tumor necrotic region (T0) measuredfrom an image, the values of the D and ρ are obtained by finding anintersection point between: (i) a line with slope ρ/D; and (ii) anisocurve corresponding to the observed equivalent radius of T0, providedby a non-linear PIHNA model.
 20. The computer program product of claim15, wherein the aggressiveness and/or prognosis of the tumor isdetermined by ρ/D.
 21. The computer program product of claim 15, whereinthe tumor is a glioma.
 22. The computer program product of claim 15,wherein the computational model comprises one or more additional termsto predict the effects of cancer therapy, to model angiogenesis withinthe tumor, or both.
 23. The computer program product of claim 15,wherein the computer program product is incorporated into imagingequipment, image analysis equipment, or both.
 24. The computer programproduct of claim 15, wherein the computer program product isincorporated into a hospital information system.
 25. The computerprogram product of claim 15, wherein the method further comprisescalculating a concentration wave-front velocity of tumor expansion fromone or more of said images.
 26. The computer program product of claim25, wherein ρ/D is calculated from an image or images obtained at twoseparate times, time₀ and time₁, using the same imaging technique. 27.The computer program product of claim 26, wherein the concentrationwave-front velocity of tumor expansion is calculated as (i) equivalentradii of a spherical tumor of volume equal to the volumes of the tumorat times time₀ and time₁ using a threshold tumor cell concentration; oras (ii) an average linear rate of expansion of the tumor alongdirections normal to the tumor's surface.
 28. The computer programproduct of claim 25, wherein the observed ratio ρ/D is obtained from alookup table.
 29. The computer program product of claim 28, wherein afirst and a second image are obtained (i) at two separate times by thesame imaging technique; or (ii) at approximately the same time usingdifferent imaging techniques.
 30. The computer program product of claim28, wherein the values for D and ρ are obtained using the concentrationwave-front velocity and observed ratio of ρ/D.
 31. A system for modelingtumor-growth in a patient, comprising: (a) one or more electroniccomputational processors; (b) one or more electronic memories; (c) thecomputer program product of claim 15; (d) a set of computer instructionsstored in the one or more electronic memories configured to execute onthe one or more electronic computational processors and generate acomputational model that estimates tumor-cell concentrations at one ormore positions within a tumor-containing tissue; and (e) an outputdevice for presenting results generated by the computational model to auser as a representation or prediction of the spatial extent of thetumor, the aggressiveness of the tumor, or both.
 32. The system of claim31 wherein the computational model comprises: (i) a cell-proliferationcomponent that includes a patient-specific tumor-cell migrationcoefficient representing a spatial rate of change of a tumor-cellconcentration gradient at a position within the tumor-containing tissue;and (ii) a diffusive-migration component that includes apatient-specific tumor-cell proliferation coefficient that represents anexponential tumor-cell proliferation rate.
 33. The system of claim 32,that further comprises: (iii) a set of computer instructions stored inthe one or more electronic memories configured to execute on the one ormore electronic computational processors and estimate thepatient-specific tumor-cell proliferation coefficient and thepatient-specific tumor-cell migration coefficient for a specific tumorwithin a specific tissue of a patient.
 34. The system of claim 33,wherein estimation of the patient-specific tumor-cell proliferationcoefficient and the patient-specific tumor-cell migration coefficient isaccomplished by calculating a first value for ρ/D for the tumor using afirst stored computational relationship between: (i) the ratio ρ/D; and(ii) two equivalent radii r₁ and r₂ of the tumor, calculated from (A) atleast two images, acquired by (1) the same or different methods, and (2)at approximately the same time, (B) threshold tumor-cell concentrationsfor the two images, and (C) constants related to initial and maximumtumor cell concentrations.
 35. The system of claim 34 wherein thecalculation of a first numerical value for ρ/D is accomplished by: (a)calculating an initial ratio of ρ/D from the first stored computationalrelationship; and (b) applying to the calculated initial ρ/D acorrection obtained from a lookup table.
 36. The system of claim 34,wherein the method further comprises the steps of: (a) determining awave front velocity of tumor expansion from at least two images acquiredby the same method at least two different times, time₀ and time₁; and(b) computing a second numerical value which is a product of thespecific tumor cell migration coefficient and the specific tumor cellproliferation coefficient using a second stored computationalrelationship relating the tumor cell migration coefficient to velocityof tumor expansion; and (c) determining a third numerical value which isthe specific tumor cell proliferation coefficient from the calculatedfirst and second numerical values; and (d) determining a fourthnumerical value which is the specific tumor cell migration coefficientfrom the computed first and second numerical values.
 37. The system ofclaim 34, wherein the method further comprises: (a) determining anequivalent radius r_(n) of necrotic tissue within the tumor from the atleast two images; (b) determining an intersection, on a two-dimensionalmatrix of pre-computed computational models, between (1) a line of slopeequal to (i) the ratio of the specific tumor-cell proliferationcoefficient to the specific tumor-cell migration coefficient and (2) anisocurve with a constant-equivalent-radius-of-necrotic-tissue value,produced by computational models of the two-dimensional matrix, equal tothe determined equivalent radius r_(n); and (c) calculating a pluralityof numerical values for the specific tumor-cell proliferationcoefficient and the specific tumor-cell migration coefficient fromindices of the computed intersection; wherein the two-dimensional matrixof pre-computed computational models comprises a plurality of tumor-cellproliferation coefficients and tumor-cell migration coefficients, eachpre-computed computational model indexed by a pair of numeric valuesrepresenting the tumor cell proliferation coefficients and the tumorcell migration coefficients.